###################################################################################
# R replication file for MAIN TEXT 
# "When do the Wealthy Support Redistribution?," BJPS, forthcoming 
# G. Feierherd, L. Schiumerini, and S. Stokes 
##################

# Clean objects and set environment options 

rm(list=ls(all=TRUE))
options(scipen=10, digits=2)
z_score = qnorm(1 - ((1 - .95)/2))

# Set colors for plots 

col1 <- "#000000"
col2 <- "#727272"
col3 <- "gray"

# Load packages 

library(foreign) 
require(sandwich) 
require(lmtest)  
library(plotrix) 
library(AER) 
library(xtable)
library(pastecs)

# Set working directory and load main file and functions 

setwd("~/")

dta <- read.dta("When do the Wealthy - Data.dta")

controldata <- dta[dta$control==1,]

source("When do the Wealthy - Functions.R") # Loads functions (explanation is provided in the source file.)


# Table 3: Number of Observations per Experimental Group -----------------------

dta$treatments <- ifelse(dta$control==1, "control", 
                          ifelse(dta$placebo==1, "Placebo", 
                                 ifelse(dta$injust==1, "Intra-just", 
                                        ifelse(dta$just==1, "Cross-class", NA))))

table(dta$nat, dta$treatments)

# Figure 2: Testing Hypothesis 1: Pocketbook Effects  --------------------------

pdf(file="plot_natural.pdf",family="Palatino", width=7,height=5.5)
b<-ivreg(controldata$p2_1z~controldata$lost,~controldata$nat)
b1<-b$coefficients[1]
b2<-b$coefficients[1]+b$coefficients[2]
bse<- sqrt(vcovHC(b, type = "HC2")[2,2])
bse2<-std.error(controldata$p2_1z[controldata$nat==1 & controldata$lost==1])
bse2  

c<-ivreg(controldata$p2_2z~controldata$lost,~controldata$nat)
c1<--c$coefficients[1]
c2<--c$coefficients[1]+c$coefficients[2]
cse<- sqrt(vcovHC(c, type = "HC2")[2,2]) 
cse2<-std.error(controldata$p2_2z[controldata$nat==1 & controldata$lost==1])
  
# Upper and lower confidence bounds

bu11 <- b1 + z_score*bse
bu21 <- b2 + z_score*bse2
bl11 <- b1 - z_score*bse
bl21 <-b2 - z_score*bse2
  
cu11 <- c1 + z_score*cse
cu21 <- c2 + z_score*cse2
cl11 <- c1 - z_score*cse
cl21 <- c2 - z_score*cse2
  
# Initialize plotting window  
  
plot(x=c(), y=c(), ylim=c(-.75, .75), xlim=c(0,8), xlab="", 
       ylab="Standardized mean",  xaxt="n", type = "n", bty='n', axes=T)

segments(1.05,b1,2.95,b2, lwd=4, lty=3, col=col3)
segments(5.05,c1,6.95,c2, lwd=4, lty=3, col=col3)

# Plot points of estimated effects
points(x=1, y=b1, pch=16, col=col1, cex=2)
points(x=3, y=b2, pch=15, col=col2, cex=2)
points(x=5, y=c1, pch=16, col=col1, cex=2)
points(x=7, y=c2, pch=15, col=col2, cex=2)
  
# Plot lines of confidence intervals
lines(x=c(1,1), y=c(bl11, bu11), lty=1, col=col1,lwd=3)
lines(x=c(3,3), y=c(bl21, bu21), col=col2, lty=1, lwd=3)
lines(x=c(5,5), y=c(cl11, cu11), lty=1, col=col1,lwd=3)
lines(x=c(7,7), y=c(cl21, cu21), col=col2, lty=1, lwd=3)
  
# Add a dashed horizontal line for zero
abline(h=0, lty=3)
  
axis(1,at=c(2,6),labels=c("Redistribution",
                                "Unemployment\ninsurance"), mgp=c(3, 2, 0))
# Insert text 
text(1-.5, b1 ,round(b1,digits=2),srt=0, cex=1)
text(3+.5, b2 ,round(b2,digits=2),srt=0, cex=1)
text(5-.5, c1 ,round(c1,digits=2),srt=0, cex=1)
text(7+.5, c2 ,round(c2,digits=2),srt=0, cex=1)

# Marginal effects: 
sig(ate=b$coefficients[2],se=bse,x=2,y=.75,f=1,est="LATE = ")
sig(ate=c$coefficients[2],se=cse,x=6,y=.75,f=1,est="LATE = ")

legend(x=1,y=-.4,c("No price hike","Price hike"), 
       lty=1, col=c(col1, col2), cex=1, seg.len=0.8,
       text.width=.6, pch=c(16,15), bty = "n")
dev.off()


# Figure 3: Testing Hypothesis2: Resentment ------------------------------------

pdf(file="plot_disadv.pdf",family="Palatino", width=8.5,height=5.5, colormodel=col3)

pol <- model_iv(dta,dta$p1_1z,dta$injust,dta$placebo,dta$lost,dta$nat) # Sames as: with(dta[dta$placebo==1 | dta$injust==1,], summary(ivreg(p1_1z~injust*lost,~injust*nat))) 
pu21 <- pol[3,1] + z_score*pol[4,1]
pu22 <- pol[3,2] + z_score*pol[4,2]
pl21 <- pol[3,1] - z_score*pol[4,1]
pl22 <- pol[3,2] - z_score*pol[4,2]
  
red <- model_iv(dta,dta$p2_1z,dta$injust,dta$placebo,dta$lost,dta$nat)
ru21 <- red[3,1] + z_score*red[4,1]
ru22 <- red[3,2] + z_score*red[4,2]
rl21 <- red[3,1] - z_score*red[4,1]
rl22 <- red[3,2] - z_score*red[4,2]
  
une <-  model_iv(dta,dta$p2_2z,dta$injust,dta$placebo,dta$lost,dta$nat)
uu21 <- une[3,1] + z_score*une[4,1]
uu22 <- une[3,2] + z_score*une[4,2]
ul21 <- une[3,1] - z_score*une[4,1]
ul22 <- une[3,2] - z_score*une[4,2]
  
# Initialize plotting window  

plot(x=c(), y=c(), ylim=c(-.75, .75), xlim=c(0,12), xlab="", 
     ylab="Standardized mean",  xaxt="n", type = "n", bty='n', axes=T)

abline(h=0, lty=3)

axis(1,at=c(2,6,10),labels=c("Opinion about\n the policy","Redistribution",
                             "Unemployment\ninsurance"), mgp=c(3, 2, 0))
# Price Hike 

segments(1.05,pol[3,1],2.95,pol[3,2], lwd=4, lty=3, col=col3)

lines(x=c(1,1), y=c(pl21, pu21), col=col1, lty=1, lwd=3)
lines(x=c(3,3), y=c(pl22, pu22), lty=1, col=col2, lwd=3)

points(x=1, y=pol[3,1], pch=16, col=col1, cex=2)
points(x=3, y=pol[3,2], pch=15, col=col2, cex=2)

text(1-.5, pol[3,1], round(pol[3,1],2),srt=0, cex=1)
text(3+.5, pol[3,2], round(pol[3,2],2),srt=0, cex=1)

sig(ate=pol[3,3],se=pol[4,3],x=1,y=.75,f=1,est="CATE = ")
  
# Redistribution

segments(5.05,red[3,1],6.95,red[3,2], lwd=4, lty=3, col=col3)

lines(x=c(5,5), y=c(rl21, ru21), col=col1, lty=1, lwd=3)
lines(x=c(7,7), y=c(rl22, ru22), lty=1, col=col2, lwd=3)

points(x=5, y=red[3,1], pch=16, col=col1, cex=2)
points(x=7, y=red[3,2], pch=15, col=col2, cex=2)

text(5-.5, red[3,1],round(red[3,1],2),srt=0, cex=1)
text(7+.5, red[3,2],round(red[3,2],2),srt=0, cex=1)

sig(ate=red[3,3],se=red[4,3],x=6,y=.75,f=1,est="CATE = ")

# Unemployment insurance

segments(9.05,une[3,1],10.95,une[3,2], lwd=4, lty=3, col=col3)

lines(x=c(9,9), y=c(ul21, uu21), col=col1, lty=1, lwd=3)
lines(x=c(11,11), y=c(ul22, uu22), lty=1, col=col2, lwd=3)

points(x=9, y=une[3,1], pch=16, col=col1, cex=2)
points(x=11, y=une[3,2], pch=15, col=col2, cex=2)

text(9-.5, une[3,1],round(une[3,1],2),srt=0, cex=1)
text(11+.5, une[3,2],round(une[3,2],2),srt=0, cex=1)

sig(ate=une[3,3],se=une[4,3],x=10,y=.75,f=1,est="CATE = ")

legend(x=9,y=-.4,c("Treated-Placebo","Treated-intraclass\ninequality"), 
       lty=1, col=c(col1, col2), cex=1, seg.len=0.8,
       text.width=.6, pch=c(16,15), bty = "n")

dev.off()


# Figure 4: Testing Hypothesis 3: Empathy ------------------------------------

pdf(file="plot_adv.pdf",family="Palatino", width=8.5, height=5.5, colormodel=col3)

pol <- model_iv(dta,dta$p1_1z,dta$injust,dta$placebo,dta$lost,dta$nat) 

pu21 <- pol[1,1] + z_score*pol[2,1]
pu22 <- pol[1,2] + z_score*pol[2,2]
pl21 <- pol[1,1] - z_score*pol[2,1]
pl22 <- pol[1,2] - z_score*pol[2,2]

red <- model_iv(dta,dta$p2_1z,dta$injust,dta$placebo,dta$lost,dta$nat)

ru21 <- red[1,1] + z_score*red[2,1]
ru22 <- red[1,2] + z_score*red[2,2]
rl21 <- red[1,1] - z_score*red[2,1]
rl22 <- red[1,2] - z_score*red[2,2]

une <- model_iv(dta,dta$p2_2z,dta$injust,dta$placebo,dta$lost,dta$nat)

uu21 <- une[1,1] + z_score*une[2,1]
uu22 <- une[1,2] + z_score*une[2,2]
ul21 <- une[1,1] - z_score*une[2,1]
ul22 <- une[1,2] - z_score*une[2,2]

# Initialize plotting window  

plot(x=c(), y=c(), ylim=c(-.75, .75), xlim=c(0,12), xlab="", 
     ylab="Standardized mean",  xaxt="n", type = "n", bty='n', axes=T)

abline(h=0, lty=3)

axis(1,at=c(2,6,10),labels=c("Opinion about\n the policy","Redistribution",
                             "Unemployment\ninsurance"), mgp=c(3, 2, 0))
# Price Hike 

segments(1.05,pol[1,1],2.95,pol[1,2], lwd=4, lty=3, col=col3)

points(x=1, y=pol[1,1], pch=16, col=col1, cex=2)
points(x=3, y=pol[1,2], pch=15, col=col2, cex=2)

lines(x=c(1,1), y=c(pl21, pu21), col=col1, lty=1, lwd=3)
lines(x=c(3,3), y=c(pl22, pu22), lty=1, col=col2, lwd=3)

text(1-.5, pol[1,1], round(pol[1,1],2),srt=0, cex=1)
text(3+.5, pol[1,2], round(pol[1,2],2),srt=0, cex=1)

sig(ate=pol[1,3],se=pol[2,3],x=1,y=.75,f=1,est="CATE = ")

# Redistribution

segments(5.05,red[1,1],6.95,red[1,2], lwd=4, lty=3, col=col3)

points(x=5, y=red[1,1], pch=16, col=col1, cex=2)
points(x=7, y=red[1,2], pch=15, col=col2, cex=2)

lines(x=c(5,5), y=c(rl21, ru21), col=col1, lty=1, lwd=3)
lines(x=c(7,7), y=c(rl22, ru22), lty=1, col=col2, lwd=3)

text(5-.5, red[1,1],round(red[1,1],2),srt=0, cex=1)
text(7+.5, red[1,2],round(red[1,2],2),srt=0, cex=1)

sig(ate=red[1,3],se=red[2,3],x=6,y=.75,f=1,est="CATE = ")

# Unemployment insurance

segments(9.05,une[1,1],10.95,une[1,2], lwd=4, lty=3, col=col3)

points(x=9, y=une[1,1], pch=16, col=col1, cex=2)
points(x=11, y=une[1,2], pch=15, col=col2, cex=2)

lines(x=c(9,9), y=c(ul21, uu21), col=col1, lty=1, lwd=3)
lines(x=c(11,11), y=c(ul22, uu22), lty=1, col=col2, lwd=3)

text(9-.5, une[1,1],round(une[1,1],2),srt=0, cex=1)
text(11+.5, une[1,2],round(une[1,2],2),srt=0, cex=1)

sig(ate=une[1,3],se=une[2,3],x=10,y=.75,f=1,est="CATE = ")

legend(x=9,y=-.4,c("Untreated-Placebo","Untreated-intraclass\ninequality"), 
       lty=1, col=c(col1, col2), cex=1, seg.len=0.8,
       text.width=.6, pch=c(16,15), bty = "n")

dev.off()


# Figure 5: Testing Hypothesis 4: Altruism --------------------------------

pdf(file="plot_crossclass.pdf",family="Palatino", width=8.5, height=5.5, colormodel=col3)

b <- lm(p1_1z~just, data=dta[dta$placebo==1 | dta$just==1,])
b1<-b$coefficients[1]
b2<- b1+b$coefficients[2]
bse<- sqrt(vcovHC(b, type = "HC2")[2,2])
bse2<-sqrt(vcovHC(b, type = "HC2")[1,1])
bse3<-std.error(dta[dta$just==1,]$p1_1z)

c <- lm(p2_1z~just, data=dta[dta$placebo==1 | dta$just==1,])
c1<-c$coefficients[1]
c2<-c1+c$coefficients[2]
cse<- sqrt(vcovHC(c, type = "HC2")[2,2])
cse2<-sqrt(vcovHC(c, type = "HC2")[1,1])
cse3<-std.error(dta[dta$just==1,]$p2_1z)

d <- lm(p2_2z~just, data=dta[dta$placebo==1 | dta$just==1,])
d1<-d$coefficients[1]
d2<-d1+c$coefficients[2]
dse<- sqrt(vcovHC(d, type = "HC2")[2,2])
dse2<-sqrt(vcovHC(d, type = "HC2")[1,1])
dse3<-std.error(dta[dta$just==1,]$p2_2z)

pu21 <- b1 + z_score*bse
pu22 <- b2 + z_score*bse3
pl21 <- b1 - z_score*bse
pl22 <- b2 - z_score*bse3

ru21 <- c1 + z_score*cse
ru22 <- c2 + z_score*cse3
rl21 <- c1 - z_score*cse
rl22 <- c2 - z_score*cse3

uu21 <- d1 + z_score*dse
uu22 <- d2 + z_score*dse3
ul21 <- d1 - z_score*dse
ul22 <- d2 - z_score*dse3

# Initialize plotting window  

plot(x=c(), y=c(), ylim=c(-.75, .75), xlim=c(0,12), xlab="", 
     ylab="Standardized mean",  xaxt="n", type = "n", bty='n', axes=T)

abline(h=0, lty=3)

axis(1,at=c(2,6,10),labels=c("Opinion about\n the policy","Redistribution",
                             "Unemployment\ninsurance"), mgp=c(3, 2, 0))
# Price Hike 

segments(1.05,b1,2.95,b2, lwd=4, lty=3, col=col3)

points(x=1, y=b1, pch=16, col=col1, cex=2)
points(x=3, y=b2, pch=15, col=col2, cex=2)

lines(x=c(1,1), y=c(pl21, pu21), col=col1, lty=1, lwd=3)
lines(x=c(3,3), y=c(pl22, pu22), lty=1, col=col2, lwd=3)

text(1-.5, b1, round(b1,2),srt=0, cex=1)
text(3+.5, b2, round(b2,2),srt=0, cex=1)

sig(ate=b$coefficients[2],se=bse2,x=2,y=.75,f=1,est="ATE = ")

# Redistribution

segments(5.05,c1,6.95,c2, lwd=4, lty=3, col=col3)

points(x=5, y=c1, pch=16, col=col1, cex=2)
points(x=7, y=c2, pch=15, col=col2, cex=2)

lines(x=c(5,5), y=c(rl21, ru21), col=col1, lty=1, lwd=3)
lines(x=c(7,7), y=c(rl22, ru22), lty=1, col=col2, lwd=3)

text(5-.5, c1, round(c1,2),srt=0, cex=1)
text(7+.5, c2, round(c2,2),srt=0, cex=1)

sig(ate=c$coefficients[2],se=cse2,x=6,y=.75,f=1,est="ATE = ")


# Unemployment insurance

segments(9.05,c1,10.95,d2, lwd=4, lty=3, col=col3)

points(x=9, y=d1, pch=16, col=col1, cex=2)
points(x=11, y=d2, pch=15, col=col2, cex=2)

lines(x=c(9,9), y=c(ul21, uu21), col=col1, lty=1, lwd=3)
lines(x=c(11,11), y=c(ul22, uu22), lty=1, col=col2, lwd=3)

text(9-.5, d1, round(d1,2),srt=0, cex=1)
text(11+.5, d2, round(d2,2),srt=0, cex=1)

sig(ate=d$coefficients[2],se=dse2,x=10,y=.75,f=1,est="ATE = ")


legend(x=9,y=-.4,c("Placebo","Crossclass\ninequality"), 
       lty=1, col=c(col1, col2), cex=1, seg.len=0.8,
       text.width=.6, pch=c(16,15), bty = "n")

dev.off()
